Functional analysis
Princess Rodriguez
4/11/2024
Install and load required ackages
library(clusterProfiler) #BiocManager::install("clusterProfiler")
library(org.Mm.eg.db) #BiocManager::install("org.Mm.eg.db") ## this is for mouse!
library(DOSE) #BiocManager::install("DOSE")
library(ggnewscale) #install.packages("ggnewscale")
library(pathview) #BiocManager::install("pathview")
#library(ggupset) #install.packages("ggupset")
library(enrichplot) #BiocManager::install("enrichplot")
library(knitr)
#remotes::install_github("YuLab-SMU/clusterProfiler") #library(BiocManager)
#BiocManager::valid()Functional analysis
The goal of functional analysis is provide biological insight, so it’s necessary to analyze our results in the context of our experimental hypothesis: Depletion of TCF7 results in an enrichment in genes required for myeloid cell identity. Therefore, we may expect an enrichment of processes/pathways related to myeloid/innate immune responses.
GO Ontologies
The Gene Ontology (GO) knowledgebase is the world’s largest source of information on the function of genes. GO terms are organized into three independent controlled vocabularies (ontologies) in a species-independent manner:
Biological process: refers to the biological role involving the gene or gene product, and could include “transcription”, “signal transduction”, and “apoptosis”. A biological process generally involves a chemical or physical change of the starting material or input.
Molecular function: represents the biochemical activity of the gene product, such activities could include “ligand”, “GTPase”, and “transporter”.
Cellular component: refers to the location in the cell of the gene product. Cellular components could include “nucleus”, “lysosome”, and “plasma membrane”.
Each GO term has a term name (e.g. DNA repair) and a unique term accession number (GO:0005125), and a single gene product can be associated with many GO terms.
clusterProfiler
We will be using clusterProfiler to perform over-representation analysis on GO terms associated with our list of significant genes.
Over-representation (or enrichment) analysis is a statistical method that determines whether genes from pre-defined sets (ex: those beloging to a specific GO term or KEGG pathway) are present more than would be expected (over-represented) in a subset of your data. In this case, the subset is your set of under or over expressed genes.
Clusterprofiler will take as input a significant gene list and a background gene list and performs statistical enrichment analysis using hypergeometric testing. The basic arguments allow the user to select the appropriate organism and GO ontology (BP, CC, MF) to test.
Prepping the data
resdata <- read.csv("TCF1KOvsWT_all_matrix_symbol.csv", header = T)
head(resdata)## input_id baseMean log2FoldChange lfcSE stat
## 1 ENSMUSG00000000001 3904.959106 0.07596003 0.07628721 0.99571117
## 2 ENSMUSG00000000028 66.299713 0.30166718 0.32777276 0.92035463
## 3 ENSMUSG00000000037 4.938477 0.10152331 1.20464180 0.08427676
## 4 ENSMUSG00000000056 833.776305 -0.22561867 0.11117532 -2.02939519
## 5 ENSMUSG00000000078 25622.396565 0.35991087 0.21278531 1.69142726
## 6 ENSMUSG00000000085 377.598537 -0.14503206 0.15036558 -0.96452961
## pvalue padj WT_CD8_rep1 WT_CD8_rep2 WT_CD8_rep3 TCF1_KO_CD8_rep1
## 1 0.31939050 0.5934865 3.601914e+03 4013.918639 3795.816820 4075.603038
## 2 0.35738747 0.6294404 5.894950e+01 53.379113 64.922751 78.882639
## 3 0.93283639 NA 9.991441e-01 9.531984 5.410229 4.533485
## 4 0.04241805 0.1492914 9.471886e+02 850.253010 899.180096 737.144665
## 5 0.09075523 0.2604725 1.800158e+04 25351.265763 24233.498716 32210.411110
## 6 0.33478047 0.6091535 3.866688e+02 381.279377 420.915833 324.597528
## TCF1_KO_CD8_rep2 TCF1_KO_CD8_rep3 symbol
## 1 3943.778610 3998.723034 Gnai3
## 2 56.672210 84.992064 Cdc45
## 3 5.060019 4.096003 Scml2
## 4 767.098842 801.792608 Narf
## 5 27000.260027 26937.364441 Klf6
## 6 388.609440 363.520276 Scmh1Create Gene lists
To perform the over-representation analysis, we need a list of
background genes and a list of significant genes. For our background
dataset we will use all genes in our resdata results table.
For the significant gene list we will use genes with p-adjusted values
less than 0.05 and absolute log2FC of 1.
first lets prep background gene set
# we want the log2 fold change
original_gene_list <- resdata$log2FoldChange
head(original_gene_list)## [1] 0.07596003 0.30166718 0.10152331 -0.22561867 0.35991087 -0.14503206# name the vector
names(original_gene_list) <- resdata$input_id
head(original_gene_list)## ENSMUSG00000000001 ENSMUSG00000000028 ENSMUSG00000000037 ENSMUSG00000000056
## 0.07596003 0.30166718 0.10152331 -0.22561867
## ENSMUSG00000000078 ENSMUSG00000000085
## 0.35991087 -0.14503206#original_gene_list# omit any NA values
background_gene_list<-na.omit(original_gene_list)# sort the list in decreasing order (required for clusterProfiler)
background_gene_list = sort(background_gene_list, decreasing = TRUE)head(background_gene_list)## ENSMUSG00000100426 ENSMUSG00000079481 ENSMUSG00000094840 ENSMUSG00000018925
## 11.002034 8.518602 8.198348 7.755458
## ENSMUSG00000091387 ENSMUSG00000023132
## 7.649971 7.247676now lets prep significant gene set
# Extract significant results (padj < 0.05) using subset()
sig_genes_df = subset(resdata, padj < 0.05)
# From significant results, we want to filter on log2fold change
sig_genes <- sig_genes_df$log2FoldChange
# Name the vector
names(sig_genes) <- sig_genes_df$input_id
# omit NA values
sig_genes <- na.omit(sig_genes)
# filter on min log2fold change (log2FoldChange > 1)
sig_genes <- names(sig_genes)[abs(sig_genes) > 1]head(sig_genes)## [1] "ENSMUSG00000000301" "ENSMUSG00000000303" "ENSMUSG00000000489"
## [4] "ENSMUSG00000000530" "ENSMUSG00000000631" "ENSMUSG00000000817"Creating enrichGO object
There are ontology options: [“BP”, “MF”, “CC”]. The
ontargument can accept either “BP” (Biological Process), “MF” (Molecular Function), and “CC” (Cellular Component) subontologies, or “ALL” for all three.There are keytypes options: This is the source of the annotation (gene ids). The options vary for each annotation. In the example of org.Mm.eg.db, the options are:
keytypes(org.Mm.eg.db)## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
## [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
## [11] "GENETYPE" "GO" "GOALL" "IPI" "MGI"
## [16] "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID"
## [21] "PROSITE" "REFSEQ" "SYMBOL" "UNIPROT"3.The different organisms with annotation databases available to use
with for the OrgDb argument can be found below: 
#go_enrich <- enrichGO(gene = sig genes,
#universe = background genes,
#OrgDb = organism Db,
#keyType = 'ENSEMBL',
#readable = T,
#ont = "BP",
#pvalueCutoff = 0.05,
#qvalueCutoff = 0.05)This next part will take a few minutes to run!
#saveRDS(go_enrich, file = "enrichGO_cd8.rds")
go_enrich <- readRDS("enrichGO_cd8.rds")## Output results from GO analysis to a table
cluster_summary <- data.frame(go_enrich)
write.csv(cluster_summary, "clusterProfiler_TCF7ko.csv")
head(cluster_summary)## ID Description GeneRatio
## GO:0007186 GO:0007186 G protein-coupled receptor signaling pathway 57/661
## GO:0002250 GO:0002250 adaptive immune response 57/661
## GO:0045785 GO:0045785 positive regulation of cell adhesion 55/661
## GO:0050865 GO:0050865 regulation of cell activation 67/661
## GO:0032943 GO:0032943 mononuclear cell proliferation 43/661
## GO:0032944 GO:0032944 regulation of mononuclear cell proliferation 37/661
## BgRatio pvalue p.adjust qvalue
## GO:0007186 341/12632 3.388225e-15 1.464052e-11 1.161270e-11
## GO:0002250 363/12632 5.566961e-14 1.202742e-10 9.540014e-11
## GO:0045785 349/12632 1.354755e-13 1.928068e-10 1.529322e-10
## GO:0050865 485/12632 1.784835e-13 1.928068e-10 1.529322e-10
## GO:0032943 248/12632 2.700771e-12 2.262129e-09 1.794295e-09
## GO:0032944 192/12632 3.659525e-12 2.262129e-09 1.794295e-09
## geneID
## GO:0007186 Ccl3/Tulp3/Ccl4/Plek/Adcy1/Pld2/Syk/Iqgap2/App/Rgs11/Atp2b4/Rgs16/Xcl1/Car2/Gnb4/P2ry1/Adgrl2/Adgrb2/Rgs12/Trpm1/Syp/Gpr83/Pde4a/Vipr1/Acpp/Cish/Mgll/Lpar6/Dgkh/Ccl5/Ptger2/Ccr7/Gpr34/Lrp1/Ramp3/Grk3/Gng2/Ccrl2/Gpr68/Baiap3/Cxcr5/Ccr4/Entpd1/F2r/Cxcr6/Mrgpre/Ccr2/Olfr524/Cx3cr1/Rnf157/Pde5a/Adgrg3/Adgrg5/Tcp11/Fzd10/Kctd12/Pde2a
## GO:0002250 Tbx21/Cd79a/Fcer2a/Mef2c/Sh2d1a/Cd247/Slamf6/Gzmb/Gata3/Pdcd1lg2/Lyst/Syk/Cd4/Cd74/Cd7/Ctla4/Il18rap/Il18r1/Xcl1/Cr2/Card9/Il6ra/Lef1/P2rx7/Lag3/Cd19/Swap70/Cd40lg/Crtam/Nedd4/Eomes/Cd226/Fut7/H2-Aa/Unc93b1/Prf1/H2-DMa/Ccr7/Prdm1/Slamf7/Cd160/Susd4/Stx11/Foxp3/Fgl2/Cd79b/Serpina3g/Serpinb9/Cd24a/Ccr2/Themis/Cx3cr1/Gzmm/Ifng/H2-Eb1/Tnfrsf13c/H2-Ab1
## GO:0045785 Cdh1/Pdgfb/Tnfsf14/Fbln1/Slc4a1/Csf1/Itga2/Gata3/Pdcd1lg2/Cyth3/Il12rb2/Myb/Pld2/Gcnt2/Syk/Il6st/Hes1/Cd4/Cd74/Ptprj/Nrp1/Xcl1/Il2ra/Il6ra/Vcam1/Lef1/Nr4a3/Epha1/Ret/Cd40lg/Il15/Pik3r2/Smad3/Chst2/Ccl5/Hpse/Fut7/H2-Aa/H2-DMa/Sirpa/Ccr7/Cd160/Zfhx3/Adk/Foxp3/Prkce/Cd24a/Ndnf/Ccr2/Sox12/Ifng/H2-Eb1/Zbtb16/Tnfrsf13c/H2-Ab1
## GO:0050865 Pdgfb/Myo18a/Tbx21/Hmox1/Mef2c/Tnfsf14/Slc4a1/Blk/Gata3/Pdcd1lg2/Il12rb2/Lyst/Myb/Plek/Id2/Pld2/Syk/Il6st/Dlg5/Hes1/Cd200/Cdkn1a/Cd4/Cd74/Ptprj/Ctla4/Xcl1/Il2ra/Il6ra/Vcam1/Lef1/Nr4a3/Lag3/Cd22/Tyrobp/Cd19/Cd40lg/Fgfr1/Il15/Ubash3b/Crtam/Tirap/Ldlr/Cd226/Ccl5/H2-Aa/H2-DMa/Sirpa/Ccr7/Prdm1/Slamf7/Cd160/Adk/Foxp3/Fgl2/Zc3h12d/Ctla2a/Cd24a/Ccr2/Tlr6/Sox12/Pde5a/Ifng/H2-Eb1/Zbtb16/Tnfrsf13c/H2-Ab1
## GO:0032943 Cd79a/Mef2c/Slc4a1/Blk/Csf1/Slamf6/Pdcd1lg2/Il12rb2/Syk/Il6st/Dlg5/Hes1/Cdkn1a/Cd4/Cd74/Csf1r/Ctla4/Xcl1/Cr2/Il2ra/Vcam1/Lef1/P2rx7/Cd22/Tyrobp/Cd19/Cd40lg/Il15/Crtam/Tirap/Ccl5/H2-Aa/Ccr7/Prdm1/Adk/Foxp3/Zc3h12d/Cd24a/Ccr2/Pde5a/Ifng/Tnfrsf13c/H2-Ab1
## GO:0032944 Mef2c/Slc4a1/Blk/Csf1/Pdcd1lg2/Il12rb2/Syk/Il6st/Dlg5/Hes1/Cdkn1a/Cd4/Cd74/Csf1r/Ctla4/Xcl1/Il2ra/Vcam1/Cd22/Tyrobp/Cd40lg/Il15/Crtam/Tirap/Ccl5/H2-Aa/Ccr7/Prdm1/Adk/Foxp3/Zc3h12d/Cd24a/Ccr2/Pde5a/Ifng/Tnfrsf13c/H2-Ab1
## Count
## GO:0007186 57
## GO:0002250 57
## GO:0045785 55
## GO:0050865 67
## GO:0032943 43
## GO:0032944 37Visualizing clusterProfiler results
clusterProfiler has a variety of options for viewing the
over-represented GO terms.
Bar plot is the most widely used method to visualize enriched terms. It depicts the enrichment scores (e.g. p values) and gene count or ratio as bar height and color. Users can specify the number of terms (most significant) or selected terms to display via the showCategory parameter.
barplot <- barplot(go_enrich,
showCategory = 10,
title = "GO Biological Pathways",
font.size = 8)
barplot
Dot plot is similar to bar plot with the capability to encode another score as dot size. Go ahead and convert the figure into a dotplot to display the top 10 GO terms by gene ratio (# genes related to GO term / total number of sig genes).
?dotplot## Help on topic 'dotplot' was found in the following packages:
##
## Package Library
## enrichplot /users/m/m/mmg23200/R/x86_64-pc-linux-gnu-library/4.2
## clusterProfiler /users/m/m/mmg23200/R/x86_64-pc-linux-gnu-library/4.2
## graphics /usr/local/lib/R/library
## lattice /usr/local/lib/R/library
##
##
## Using the first match ...#from enrichplot librarydot <- dotplot(go_enrich,
showCategory = 10,
font.size = 8)
dot
We can also visualize enriched GO terms as a directed acyclic graph:
pathway <- goplot(go_enrich, showCategory = 10)
pathway
Gene Concept Network
Both the barplot() and dotplot() only
displayed most significant terms, but some users may want to know which
genes are involved in these significant terms.
In order to consider the potentially biological complexities in which
a gene may belong to multiple categories, we can use the
cnetplot() function. The cnetplot() depicts
the linkages of genes and biological concepts (e.g. GO terms or KEGG
pathways) as a network.
cnetplot(go_enrich,
showCategory = 5,
categorySize="pvalue",
foldChange=background_gene_list,
vertex.label.font=6,
cex_label_gene = 0.5,
cex_label_category= 0.8)
cnetplot(go_enrich,
showCategory = 5,
categorySize="pvalue",
foldChange=background_gene_list,
vertex.label.font=6,
cex_label_category = 0.5,
node_label="category")
Select which labels to be displayed. Options include ‘category’, ‘gene’, ‘all’(the default) and ‘none’.
library(ggplot2)
plot_color <- cnetplot(go_enrich,
showCategory = 5,
categorySize="pvalue",
foldChange=background_gene_list,
vertex.label.font=6,
cex_label_gene = 0.5,
cex_label_category= 0.8) +
scale_color_gradient2(name='fold change', low='darkgreen', high='firebrick')
plot_color
cnetplot(go_enrich,
showCategory = 5,
categorySize="pvalue",
foldChange= background_gene_list,
circular = TRUE,
colorEdge = TRUE,
cex_label_gene = 0.35,
vertex.label.font=6)
Circular plots can be quite powerful, if the list was smaller, here things get a bit crowded - overall makes it hard to interpret. Moral of the story, just because you can make it doesn’t mean you should!
tree_cd8 <- pairwise_termsim(go_enrich)
treeplot(tree_cd8)
KEGG Pathway Enrichment
For KEGG pathway enrichment using the gseKEGG() function, we need to convert id types. We can use the bitr function for this (included in clusterProfiler). It is normal for this call to produce some messages / warnings.
# Convert gene IDs for enrichKEGG function
# We will lose some genes here because not all IDs will be converted
ids<-bitr(names(original_gene_list),
fromType = "ENSEMBL",
toType = "ENTREZID",
OrgDb="org.Mm.eg.db")
# remove duplicate IDS (here I use "ENSEMBL", but it should be whatever was selected as keyType)
dedup_ids = ids[!duplicated(ids[c("ENSEMBL")]),]# Create a new dataframe df2 which has only the genes which were successfully mapped using the bitr function above
df2 = resdata[resdata$input_id %in% dedup_ids$ENSEMBL,]# Create a new column in df2 with the corresponding ENTREZ IDs
df2 <- merge(resdata, dedup_ids, by.x = "input_id", by.y = "ENSEMBL")
str(df2)## 'data.frame': 14850 obs. of 15 variables:
## $ input_id : chr "ENSMUSG00000000001" "ENSMUSG00000000028" "ENSMUSG00000000037" "ENSMUSG00000000056" ...
## $ baseMean : num 3904.96 66.3 4.94 833.78 25622.4 ...
## $ log2FoldChange : num 0.076 0.302 0.102 -0.226 0.36 ...
## $ lfcSE : num 0.0763 0.3278 1.2046 0.1112 0.2128 ...
## $ stat : num 0.9957 0.9204 0.0843 -2.0294 1.6914 ...
## $ pvalue : num 0.3194 0.3574 0.9328 0.0424 0.0908 ...
## $ padj : num 0.593 0.629 NA 0.149 0.26 ...
## $ WT_CD8_rep1 : num 3.60e+03 5.89e+01 9.99e-01 9.47e+02 1.80e+04 ...
## $ WT_CD8_rep2 : num 4013.92 53.38 9.53 850.25 25351.27 ...
## $ WT_CD8_rep3 : num 3795.82 64.92 5.41 899.18 24233.5 ...
## $ TCF1_KO_CD8_rep1: num 4075.6 78.88 4.53 737.14 32210.41 ...
## $ TCF1_KO_CD8_rep2: num 3943.78 56.67 5.06 767.1 27000.26 ...
## $ TCF1_KO_CD8_rep3: num 3998.7 85 4.1 801.8 26937.4 ...
## $ symbol : chr "Gnai3" "Cdc45" "Scml2" "Narf" ...
## $ ENTREZID : chr "14679" "12544" "107815" "67608" ...# Create a vector of the gene universe
kegg_gene_list <- df2$log2FoldChange# Name vector with ENTREZ ids
names(kegg_gene_list) <- df2$ENTREZIDhead(kegg_gene_list)## 14679 12544 107815 67608 23849 29871
## 0.07596003 0.30166718 0.10152331 -0.22561867 0.35991087 -0.14503206# omit any NA values
kegg_gene_list<-na.omit(kegg_gene_list)# sort the list in decreasing order (required for clusterProfiler)
kegg_gene_list = sort(kegg_gene_list, decreasing = TRUE)We just created the background list, similar to what we did above. Now, we will move on to create the significant list.
# Extract significant results from df2
kegg_sig_genes_df = subset(df2, padj < 0.05)# From significant results, we want to filter on log2fold change
kegg_sig_genes <- kegg_sig_genes_df$log2FoldChange# Name the vector with the ENTREZID
names(kegg_sig_genes) <- kegg_sig_genes_df$ENTREZIDhead(kegg_sig_genes)## 74204 14673 18618 12550 67027 54204
## -0.3414865 -0.3829829 -1.2653917 4.2399463 0.3273819 -0.2018348# omit NA values
kegg_sig_genes <- na.omit(kegg_sig_genes)# filter on log2fold change
kegg_sig_genes <- names(kegg_sig_genes)[abs(kegg_sig_genes) > 1]BTW, the clusterProfiler package provides
search_kegg_organism() function to help searching supported
organisms.
sal_ent <- search_kegg_organism('Salmonella enterica', by='scientific_name')
dim(sal_ent)## [1] 50 3#sal_entCreating an enrichKEGG object
organism KEGG Organism Code: The full list is here: https://www.genome.jp/kegg/catalog/org_list.html (need the 3 letter code). I define this as kegg_organism first, because it is used again below when making the pathview plots.
keyType one of ‘kegg’, ‘ncbi-geneid’, ‘ncib-proteinid’ or ‘uniprot’.
head(kk)## ID
## mmu04060 mmu04060
## mmu04061 mmu04061
## mmu04640 mmu04640
## mmu05321 mmu05321
## mmu05323 mmu05323
## mmu05200 mmu05200
## Description
## mmu04060 Cytokine-cytokine receptor interaction - Mus musculus (house mouse)
## mmu04061 Viral protein interaction with cytokine and cytokine receptor - Mus musculus (house mouse)
## mmu04640 Hematopoietic cell lineage - Mus musculus (house mouse)
## mmu05321 Inflammatory bowel disease - Mus musculus (house mouse)
## mmu05323 Rheumatoid arthritis - Mus musculus (house mouse)
## mmu05200 Pathways in cancer - Mus musculus (house mouse)
## GeneRatio BgRatio pvalue p.adjust qvalue
## mmu04060 30/320 129/5552 1.965334e-11 4.451543e-09 3.867418e-09
## mmu04061 18/320 47/5552 3.049002e-11 4.451543e-09 3.867418e-09
## mmu04640 18/320 67/5552 2.161501e-08 2.103861e-06 1.827796e-06
## mmu05321 14/320 46/5552 1.502688e-07 1.096962e-05 9.530204e-06
## mmu05323 12/320 56/5552 6.128943e-05 3.181992e-03 2.764456e-03
## mmu05200 40/320 372/5552 7.111370e-05 3.181992e-03 2.764456e-03
## geneID
## mmu04060 11482/14103/20302/50930/12977/16162/20303/16195/12504/12978/16174/16182/16963/16184/16194/21949/21947/16168/16154/20304/12775/12145/100532/12773/80901/12772/13051/15978/72049/16188
## mmu04061 20302/50930/12977/20303/16195/12978/16174/16182/16963/16184/16194/16154/20304/12775/12145/12773/12772/13051
## mmu04640 14128/12977/16398/12504/12978/12482/12516/12902/16184/16194/12483/12478/14960/14998/12484/14969/16188/14961
## mmu05321 57765/14462/16162/16174/16182/24088/17127/14960/14998/20371/15978/17132/14969/14961
## mmu05323 20302/12977/11975/12477/24088/16168/20304/14960/14998/15978/14969/14961
## mmu05200 12550/18591/14103/15368/16001/17873/16398/16162/432530/18806/69635/16195/15205/12575/18712/58801/12978/21416/226519/16184/13555/14696/16194/16842/19713/14182/16168/18709/13143/17127/67168/19217/18131/14702/14062/15978/21415/235320/16188/93897
## Count
## mmu04060 30
## mmu04061 18
## mmu04640 18
## mmu05321 14
## mmu05323 12
## mmu05200 40saveRDS(kk, file = "enrichKEGG_cd8.rds")
kk <- readRDS("enrichKEGG_cd8.rds")#saveRDS(gseaKEGG_results, file = "gseaKEGG_results.rds")
gseaKEGG_results <- readRDS("gseaKEGG_results.rds")Visualize enriched KEGG pathways
To view the KEGG pathway, users can use the browseKEGG function, which will open a web browser and highlight enriched genes.
#browseKEGG(kk, 'mmu05202')library("pathview")
mmu05152 <- pathview(gene.data = kegg_gene_list,
pathway.id = "mmu05152",
species = "mmu",
limit = list(gene=max(abs(kegg_gene_list)), cpd=1))
knitr::include_graphics("mmu05152.pathview.png")
Session info
sessionInfo()## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggplot2_3.5.0 knitr_1.46 enrichplot_1.18.4
## [4] pathview_1.38.0 ggnewscale_0.4.10 DOSE_3.24.2
## [7] org.Mm.eg.db_3.16.0 AnnotationDbi_1.60.2 IRanges_2.32.0
## [10] S4Vectors_0.36.2 Biobase_2.58.0 BiocGenerics_0.44.0
## [13] clusterProfiler_4.6.2
##
## loaded via a namespace (and not attached):
## [1] fgsea_1.24.0 colorspace_2.1-0 ggtree_3.6.2
## [4] gson_0.1.0 qvalue_2.30.0 XVector_0.38.0
## [7] fs_1.6.3 aplot_0.2.2 rstudioapi_0.16.0
## [10] farver_2.1.1 graphlayouts_1.1.1 ggrepel_0.9.5
## [13] bit64_4.0.5 fansi_1.0.6 scatterpie_0.2.2
## [16] codetools_0.2-18 splines_4.2.2 cachem_1.0.8
## [19] GOSemSim_2.24.0 polyclip_1.10-6 jsonlite_1.8.8
## [22] GO.db_3.16.0 png_0.1-8 graph_1.76.0
## [25] ggforce_0.4.2 compiler_4.2.2 httr_1.4.7
## [28] Matrix_1.5-1 fastmap_1.1.1 lazyeval_0.2.2
## [31] cli_3.6.2 tweenr_2.0.3 htmltools_0.5.8.1
## [34] tools_4.2.2 igraph_2.0.3 gtable_0.3.4
## [37] glue_1.7.0 GenomeInfoDbData_1.2.9 reshape2_1.4.4
## [40] dplyr_1.1.4 fastmatch_1.1-4 Rcpp_1.0.12
## [43] jquerylib_0.1.4 vctrs_0.6.5 Biostrings_2.66.0
## [46] ape_5.7-1 nlme_3.1-160 ggraph_2.2.1
## [49] xfun_0.43 stringr_1.5.1 lifecycle_1.0.4
## [52] XML_3.99-0.16.1 org.Hs.eg.db_3.16.0 zlibbioc_1.44.0
## [55] MASS_7.3-58.1 scales_1.3.0 tidygraph_1.3.1
## [58] parallel_4.2.2 KEGGgraph_1.58.3 RColorBrewer_1.1-3
## [61] yaml_2.3.8 memoise_2.0.1 gridExtra_2.3
## [64] downloader_0.4 ggfun_0.1.4 HDO.db_0.99.1
## [67] yulab.utils_0.1.4 sass_0.4.9 stringi_1.8.3
## [70] RSQLite_2.3.6 highr_0.10 tidytree_0.4.6
## [73] BiocParallel_1.32.6 GenomeInfoDb_1.34.9 rlang_1.1.3
## [76] pkgconfig_2.0.3 bitops_1.0-7 evaluate_0.23
## [79] lattice_0.20-45 purrr_1.0.2 labeling_0.4.3
## [82] treeio_1.22.0 patchwork_1.2.0 cowplot_1.1.3
## [85] shadowtext_0.1.3 bit_4.0.5 tidyselect_1.2.1
## [88] plyr_1.8.9 magrittr_2.0.3 R6_2.5.1
## [91] generics_0.1.3 DBI_1.2.2 pillar_1.9.0
## [94] withr_3.0.0 KEGGREST_1.38.0 RCurl_1.98-1.14
## [97] tibble_3.2.1 crayon_1.5.2 utf8_1.2.4
## [100] rmarkdown_2.26 viridis_0.6.5 grid_4.2.2
## [103] data.table_1.15.4 Rgraphviz_2.42.0 blob_1.2.4
## [106] digest_0.6.35 tidyr_1.3.1 gridGraphics_0.5-1
## [109] munsell_0.5.1 viridisLite_0.4.2 ggplotify_0.1.2
## [112] bslib_0.7.0